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ABSTRACT 



Precise radial-velocity measurements for data acquired with the HARPS spectrograph infer that three planets orbit the M4 dwarf 
star GJ876. In particular, we confirm the existence of planet d, which orbits every 1.93785 days. We find that its orbit may have 
significant eccentricity (e=0.14), and deduce a more accurate estimate of its minimum mass of 6.3 M s . Dynamical modeling of the 
HARPS measurements combined with literature velocities from the Keck Observatory strongly constrain the orbital inclinations of 
the b and c planets. We find that i h = 48.9° ± 1.0° and i c = 48. 1° ± 2. 1°, which infers the true planet masses of M b = 2.64 ± 0.04 M ]up 
and M c = 0.83 ± 0.03 M lup , respectively. Radial velocities alone, in this favorable case, can therefore fully determine the orbital 
architecture of a multi-planet system, without the input from astrometry or transits. 

The orbits of the two giant planets are nearly coplanar, and their 2: 1 mean motion resonance ensures stability over at least 5 Gyr. The 
libration amplitude is smaller than 2°, suggesting that it was damped by some dissipative process during planet formation. The system 
has space for a stable fourth planet in a 4:1 mean motion resonance with planet b, with a period around 15 days. The radial velocity 
measurements constrain the mass of this possible additional planet to be at most that of the Earth. 

Key words, stars: individual: GJ 876 - stars: planetary systems - techniques: radial velocities - methods: observational 



1. Introduction 

Also known as IL Aqr, GJ 876 is only 4.72 pc away from our 
Sun and is thus the closest star known to harbor a multi-planet 
system. It has been tracked by several instruments and tele- 
scopes from the very beginning of the planetary hunt in 1994, 
namely the Hamilton (Lick), HIRES (Keck) echelle spectrome- 
ters, ELODIE (Haute-Provence), and CORALIE (La Silla) spec- 
trographs. More recently, it was also included in the HARPS 
program. 

The HARPS search for southern extra-solar planets is an 
extensive radial-velocity survey of some 2000 stars in the so- 
lar neighborhood conducted with the HARPS spectrograph on 
the ESO 3.6-m telescope at La Silla (Chile) in the framework 
of Guaranteed Tim e Observations granted to the HARPS build- 
ing consortium dMavor et alJl2003l> . About 10% of the HARPS 
GTO time were dedicated to observe a volume-limited sample of 
~1 10 M dwarfs. T his program has proven to be very efficient in 
finding Neptunes (iBonfils et alJl20 05b. 2007) and Super-Earths 



* Based on observations made with the HARPS instrument on 
the ESO 3.6 m telescope at La Silla Observatory under the GTO 
programme ID 072.C-0488; and on observations obtained at the 
Keck Observatory, which is operated jointly by the University of 
California and the California Institute of Technology. The table with 
the HARPS radial velocities is available in electronic form at the 
CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via 



http://cdsweb.u-strasbg.fr/cgi-bin/qcat?J/A+A/?? 
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(lUdrv et alJ2007tlForveille et al.l2009r.lMavor et al.ll2 009) down 
to m sin i = 1.9 M m . Because M dwarfs are more favorable tar- 
gets to searches for lower mass and/or cooler planets than around 
Sun-like stars, the observational effort dedicated by ESO has in- 
creased by a factor of four, to an extended sample of 300 stars. 

The first planet around GJ 876, a Jupiter-mass planet 
with a period of about 61 day s, was simultaneo us re- 
ported bvlDelfosse et al.1 d!998l) and lMarcv etal] d!998l) . Later, 
Mar cv et alJ (1200 ll) found that the 61 -day signal was produced 
by two planets in a 2:1 mean motion resonance, the inner 
one with 30-day period also being a gas giant. While inves- 
tigating the dynam ical interactions between those two planets, 
iRivera et al.1 0005) found evidence of a third planet, with an or- 
bital period less than 2 days and a minimum mass of 5.9 M e , 
which at that time was the lowest mass detected companion to a 
main-sequence star other than the Sun. 

Systems with two or more interacting planets dramatically 
improve our ability to understand planetary formation and evolu- 
tion, since dynamical analysis can both constrain their evolution- 
ary history and more accurately determine their orbital "struc- 
ture". Amongst known multi-planet systems, planets GJ876/? 
and c have by far the strongest mutual gravitational interactions, 
and all their orbital quantities change quite rapidly. Since the 
radial velocity of GJ876 has been monitored for the past 15 
years, the true masses of these two planets can be determined 
by adjusting the inclination of their orbital planes when fitting 
dynamical orbits to the observational data. This was first at- 
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tempted by iRivera & Lissauerl (1200 lh . who found a broad min- 
imum in the residuals to the observed radial velocities for in- 
clinations higher than 30°, the best-fit function corresponding 
to ~ 37°. Astrometric observations of GJ876 with the Hubble 
Space Telescope subs e quently suggested a higher value of ~ 85° 
dBenedict et al.ll2002l) . IRivera et al] d2005l) re-examined the dy- 
namical interactions with many additional Keck radial velocity 
measurements, and fou nd an intermediate inclination of ~ 50°. 
Bean & Seifahrt (2009) reconciled the astrometry and radial ve- 
locities by performing a joint adjustment to the Keck and HST 
datasets, and showed that both are consistent with ~ 50°. Early 
inclination determinations were, in retrospect, affected by small- 
number statistics (for astrometry) and by a modest signal-to- 
noise ratio in the radial velocity residuals. 

In the present study, we reanalyze the GJ 876 system, includ- 
ing 52 additional high precision radial velocity measurements 
taken with the HARPS spectrograph. We confirm the presence of 
a small planet d and determine the masses of b and c with great 
precision. Section 2 summarizes information about the GJ876 
star. Its derived orbital solution is described in Sect. 3, and the 
dynamical analysis of the system is discussed in Sect. 4. Finally, 
conclusions are drawn in Sect. 5. 

2. Stellar characteristics of GJ876 

GJ876 (ILAq uarii, Ross 780, HIP 113020) is a M4 dwarf 
dHawlev et alj 119961) in the Aquari us constellation. At 4.7 pc 
(tt = 2 1 2.69 + 2. 1 mas - lESAl 1997b . it is the 4 1 st closest known 
stellar systerrQ and only the 3 rd closest known planetary system 
(after e Eridani, and slightly further away than GJ 674). 

Its photometry (V = 10.162 + 0.009, K = 5.010 + 0.021 - 
iTuron et aDll993HCutri et al.ll2003l) and its parallax imply abso- 
lute magnitudes of M v - 1 1.80 ± 0.0 5 and M K = 6.65 + 0.05. 
The J - K co l or of GJ 876 (= 0.92 - ICutri et al.ll2003l) and the 
Legg ett et al.l (1200 lh color-bolometric relation infer a K-band 
bolometric correction of BCk = 2.80, and a .013 L n lumi- 
nosity . The K-band mass-luminosity relation of iDelfosse et alj 
(2000) implies a 0.334M o mass, wh ich is comparable to the 
0.32Mm derived bvlRiyera et alJ d2005h from lHenrv & McCarthy! 
(fl993b . iBonfils et all (2005a) estimate it has an approximately 
solar metallicity ([Fe/H] = 0.02 ± 0.21. lJohnson & Ad ds (2009) 
revised this metallicity calibration for the most metal-rich M 
dwarfs and found that GJ876 has an above solar metallicity 
([Fe/H] > +0.3dex). 

In terms of magnetic activity, GJ 876 is also an average star 
in our sample regarding its . Its Ca n H&K emission is almost 
twice the emission of Barnard's star, an old star in our sample of 
the same spectral type, but still comparable to many stars with 
low jitter (rms~ 2 ms~'). On the one hand, its long rotational 
period and magnetic activity imply an old age (> 0.1 Gyr). On 
the other hand, its UV W Galactic velo cities place GJ 876 in the 
young disk population (LeggetjQ992), suggesting an age < 5 
Gy r, and therefore brac keting its age to ~ 0.1 - 5 Gyr. 

IRivera et al.l d2005i) monitored GJ 876 for photometric vari- 
ability and found a 96.7-day periodic variation with a ~1% am- 
plitude, hence identifying the stellar rotation. This short rota- 
tional period is particularly helpful for radial-velocity measure- 
ments as a dark spot covering 1% of GJ 876's surface and lo- 
cated on its equator would produce a Doppler modulation with a 
maximum semi-amplitude of only ~1.5 m/s. 

Finally, t he high proper motion of GJ 876 (1.17 arcsecyr 1 
- lESAll 19971) changes the orientation of its velocity vector along 



Table 1. Observed and inferred stellar parameters of GJ 876. 



Parameter 




GJ876 


Sp 




M4V 


V 


[mag] 


10.162 


n 


[mas] 


212.69 ± 2.10 


M v 


[mag] 


11.80 ±0.05 


[Fe/H] 


[dex] 


0.05 ± 0.20 


L 


[L a ] 


0.013 


M, 


[M Q ] 


0.334 ± 10% 


Plot 


[day] 


96.7 


v sin i 


[km s -1 ] 


0.16 


age(log^ K ) 


[Gyr] 


0.1-5 



http://www.chara.gsu.edu/RECONS/TOP100.posted.htm 



the line of sight (e.g. lKiirster et al.l2003l) resulting in an apparent 
secular acceleration of 0.15 ms^'yr -1 . Before our orbital analy- 
sis, we removed this drift for both HARPS and already published 
data. 



3. Orbital solution for the GJ876 system 

3.1. Keck + HARPS 

The HARPS observations of GJ 876 started in December 2003 
and have continued for more than four years. During this time, 
we acquired 52 radial-velocity measurements with an average 
precision of 1.0 m/s. 

Before the HARPS program, the star GJ 876 had already 
been followed since June 1997 by the HIRES spectrometer 
mounted on the 10-m telescope I at Keck Observatory (Hawaii). 
The last published set of data acquired at Keck are from 
December 2004 and correspond to 155 radi al-velocity measure - 
ments with an average precision of 4.3 m/s (IRivera et al.l l2005). 

When combining the Keck and the HARPS data for GJ 876, 
the time span for the observations increases to more than 11 
years, and the secular dynamics of the system can be far more 
tightly constrained, although the average precision of the Keck 
measurements is about four times less accurate than for HARPS. 



3.2. Two planet solution 

With 207 measurements (155 from Keck and 52 from HARPS), 
we are now able to determine the nature of the three bodies in 
the system with gre at accuracy. Using the iterative Levenberg- 
Marquardt method dPress et al.l [T992), we first attempt to fit 
the complete set of radial velocities with a 3 -body Newtonian 
model (two planets) assuming coplanar motion perpendicular to 
the plane o f the sky, similarly to that achieved for the system 
HD 202206 dCorreia et al.ll2005h . This fit yields the well-known 
planetary companions at P/, = 61 day, and P c = 30 day, and an 
adjustment of yjx*- = 2.64 and rms = 4.52m/s. 

3.3. Fitting the inclinations 

Because of the proximity of the two planets and their high min- 
imum masses, the gravitational interactions between these two 
bodies are strong. Since the observational data cover more than 
one decade, we can expect to constrain the inclination of their or- 
bital planes. Without loss of generality, we use the plane of the 
sky as a reference plane, and choose the line of nodes of planet 
b as a reference direction in this plane (Q/, = 0). We thus add 
only three free parameters to our fit: the longitude of the node 
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KECK + HARPS 




50500 51000 51500 52000 52500 53000 53500 54000 54500 
JD-2400000 [days] 

Fig. 1. Keck and HARPS radial velocities for GJ 876, superim- 
posed on the 3-body Newtonian (two planets) orbital solution for 
planets b and c (orbital parameters taken from Tabled . 



Q. c of planet c, and the inclinations if, and i c of planets b and c 
with respect to the plane of the sky. 

This fit yields if, = 51°, i c = 63°, and Q c = -1°, and an ad- 
justment of -\/x^ - 2.57 and rms - 4.35m/s. Although there is 
no significant improvement to the fit, an important difference ex- 
ists: the new orbital parameters for both planets show some devi- 
ations and, more importantly, the inclinations of both planets are 
well constrained. In Fig.Q] we plot the Keck and HARPS radial 
velocities for GJ 876, superimposed on the 3-body Newtonian 
orbital solution for planets b and c. 

3.4. Third planet solution 

The residuals around the best-fit two-planet solution are small, 
but still larger than the internal errors (Fig.Q]i. We may then 
search for other companions with different orbital periods. 
Performing a frequency analysis of the velocity residuals (FigJTJ, 
we find an important peak signature at P — 1.9379 day (Fig|2]>. 
This peak was already presen t in the Keck data analysis per- 
formed by iRivera et al.l d2005l) . but it is here reinforced by the 
HARPS data. In Fig. [2] we also plot the window function and 
conclude that the above mention peak cannot be an aliasing of 
the observational data. 

To test the planetary nature of the signal, we performed a 
Keplerian fit to the residuals of the two planets. We found an 
elliptical orbit with P d = 1.9378 day, e d - 0.14, and M d sin i d = 
5.4 M e . The adjustment gives y/x* = 1.52 and rms = 2.63 m/s, 
which represents a substantial improvement with respect to the 
system with only two companions (Fig.[T|i. 

3.5. Complete orbital solution 

Starting with the orbital parameters derived in the above sec- 
tions, we can now consider the best-fit orbital solution for 
GJ 876. Using the iterative Levenberg-Marquardt method, we 
then fit the complete set of radial velocities with a 4-body 
Newtonian model. In this model, 20 parameters out of 23 possi- 
ble are free to vary. The three fixed parameters are Qb = 0° (by 




10 2 
Period [day] 

Fig. 2. Frequency analysis and window function for the two- 
planet Keck and HARPS residual radial velocities of GJ 876 
(FigHJ. An important peak is detected at P = 1 .9379 day, which 
can be interpreted as a third planetary companion in the system. 
Looking at the window function, we can see that this peak is not 
an artifact. 



definition of the reference plane), but also Cl d = 0° and i d = 50°, 
because the current precision and time span of the observations 
are not large enough to constrain these two orbital parameters. 
However, the quality of the fit does not change if we choose other 
values for these parameters. 

The orbital parameters corresponding to the best-fit solution 
are listed in Table|2] In particular, we find if, = 48.9° + 1.0° and 
i c =48.1° ±2.1°, which infer, respectively, M b = 2.64±0.04M Jup 
andM f = 0.83+0.03 Mj up for the true masses of the planets. This 

fit yields an adjustment of -\[x^ = 1-37 and rms = 2.29 ms _1 , 
which represents a significant improvement with respect to all 
previous solutions. We note that the predominant uncertainties 
are related to the star's mass (Table[TJ, but these are not folded 
into the quoted error bars. 

The best-fit orbit for planet d is also eccentric, in contrast 
to previous determinations in which its value was c onstrained 
to be zero dRivera et al J 120051; iBean & Seifahrtll2009l) . We can 
also fix this parameter and our revised solution corresponds to an 
adjustment with yfx* = 1-40 and rws = 2.34m/s. These values 
are compatible with the best-fit solution in Table|2] so we cannot 
rule out the possibility that future observations will decrease the 
eccentricity of planet d to a value close to zero. 
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Table 2. Orbital parameters for the planets orbiting GJ 876, obtained with a 4-body Newtonian fit to observational data from Keck 
and HARPS. 



Param. 


[unit] 


b 


c 


d 


Date 


[JD] 




2455 000.00 (fixed) 




V(Keck) 


[km/s] 




0.0130 ± 0.0004 




V (HARPS) 


[KITl/S J 




-1.3388 ±0.0004 




P 


[day] 


61.067 + 0.011 


30.258 + 0.009 1.93785 ± 0.00002 


A 


[deg] 


35.61 +0.14 


158.62 ±0.80 


29.94 ± 3.30 


e 




0.029 ± 0.001 


0.266 ± 0.003 


0.139 ±0.032 


OJ 


[deg] 


275.52 ± 2.67 


275.26 ± 1.25 


170.60 ± 15.52 


A 


Lm/sJ 


Z1Z.Z4 ± U.jj 


86.15 ±0.40 


6.67 ± 0.26 




t uc gJ 


^ro.y~j m \j.y J 


48.07 ± 2.06 


50 (fixed) 


Q. 


[deel 


(fixed) 


-2.32 ± 0.94 


(fixed) 


a t sin i 


[10 3 AU] 


1.19 


0.23 


1.2 x 10~ 3 


f(M) 


[10 M Q ] 


60.41 


1.79 


5.8 x 10" 5 


M sin ! 


[M ffi ] 






6.3 


M 


[M Jup ] 


2.64 


0.83 




a 


[AU] 


0.211 


0.132 


0.021 


^meas 






207 




Span 


[day] 




4103 




V? 






1.37 




rms (Keck) 


[m/s] 




4.25 




rms {HARPS) 


[m/s] 




1.80 





Errors are given by the standard deviation tr, and A is the mean longitude of the date. 



are not as strong as the mutual interactions between the two out- 
ermost planets. 

We can nevertheless estimate the time span of GJ 876 radial 
velocity measures that is necessary before we will be able to de- 
termine the inclination id with good accuracy. For that purpose, 
we fit the observational data for different fixed values of the in- 
clination of planet d (i d = 10°, 20°, 50°, and 90°). In Fig.g] we 
plot the evolution of the radial velocity differences between each 
solution and the solution at 90°. We also plot the velocity resid- 
uals of the best-fit solution from Table[2] We observe that, with 
the current HARPS precision, the inclination will be constrained 
around 2025 if id is close to 10°, but if this planet is also copla- 
nar with the other two (id ~ 50°), its precise value cannot be 
obtained before 2050 (Fig.©. 

3.7. Other instruments 

Besides the Keck and HARPS programs, the star GJ 876 was 
also followed by many other instruments. The oldest observa- 
tional data were acquired using the Hamilton echelle spectrome- 
ter mounted on the 3-m Shane telescope at the Lick Observatory. 
The star was followed from November 1994 until December 
2000 and 16 radial velocity measurements were acquired with 
an average precision of 25 m/s dMarcv etalj2 001). During 1998, 
between July and September a quick series of 40 radial velocity 
observations were performed using the CORALIE echelle spec- 
trograph mounted on the 1.2-m S wiss telescope at La Silla, with 
an average precision of 30 m/s dDelfosse et alJll998l) . Finally, 
from October 1995 to October 2003, the star was also observed 
at the Haute-Provence Observatory (OHP, France) using the 
ELODIE high-precision fiber-fed echelle spectrograph mounted 
at the Cassegrain focus of the 1.93-m telescope. Forty-six radial 
velocity measurements were taken with an average precision of 
18 m/s (data also provided with this paper via CDS). 

We did not consider these data sets in the previous anal- 
ysis because their inclusion would have been more distracting 



Gl 876 d KECK + HARPS 




0.2 0.4 0.6 0.8 1 



20 - 




50500 51000 51500 52000 52500 53000 53500 54000 54500 



JD-2400000 [days] 

Fig. 3. Phase-folded residual radial velocities for GJ 876 when 
the contributions from the planets b and c are subtracted. Data 
are superimposed on a Keplerian orbit of P = 1 .9378 day and 
e = 0.14 (the complete set of orbital parameters are those from 
Table|2]i. The respective residuals as a function of Julian date are 
displayed in the lower panel. 



3.6. Inclination of planet d 

With the presently available data, we were able to obtain with 
great accuracy the inclinations of planets b and c. However, this 
was not the case for planet d, whose inclination was held fixed 
at 50° in the best-fit solution (Tabled. We are still unable to 
determine the inclination of the innermost planet, because the 
gravitational interactions between this planet and the other two 
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Table 3. Orbital parameters for the planets orbiting GJ 876, obtained with a 4-body Newtonian fit to all available radial velocity 
observational data (Lick, Keck, CORALIE, ELODIE, and HARPS). 



Param. 


[unit] 


b 


c 


d 


Date 


[JD] 




2455 000.00 (fixed) 






[km/s] 




-0.0304 ± 0.0058 




VKeck 


[km/s] 




0.0131 ±0.0058 




^CORALIE 


[km/s] 




-1.8990 ±0.0061 




V ' ELODIE 


[km/s] 




-1.8624 ±0.0064 




Vharps 


[km/s] 




-1.3389 ±0.0058 




P 


[day] 


61.065 ±0.012 


30.259 ±0.010 1 .93785 ± 0.00002 


A 


[dog] 


35.56 ±0.15 


159.05 ±0.76 


29.37 ±3.21 


e 




0.031 ±0.001 


0.265 ± 0.002 


0.124 ±0.032 


CO 


[deg] 


274.69 ± 2.57 


274.64 ± 1.19 


168.45 ± 17.21 


K 


[m/s] 


212.09 ± 0.33 


85.91 ± 0.40 


6.60 ± 0.26 


i 


Ldeg] 


a q no i r\ c\a 


45.67 ± 1.81 


50 (fixed) 


o 


L ae gJ 




-1.67 ±0.84 


(fixed) 


a i sin; 


[10~ 3 AU] 


1.19 


0.23 


1.2 x 10~ 3 


KM) 


[10- 9 M ] 


60.27 


1.78 


5.6 x 10~ 5 


M sin i 


[Me] 






6.2 


M 


[M Jup ] 


2.64 


0.86 




a 


[AU] 


0.211 


0.132 


0.021 








309 




Span 


[day] 




5025 




V? 






1.46 




rms 


[m/s] 




3.01 





Errors are given by the standard deviation cr, and A is the mean longitude of the date. 




-20 1 — L-t — 1 1 1 1 1 1 1 1 LJ 1 

2000 2010 2020 2030 2040 2050 

year 



Fig. 4. Radial velocity differences between orbital solutions with 
id = 10°, 20°, and 50° and the orbital solution with i d = 90°. 
We also plot the velocity residuals of the best-fit solution from 
Table[2] We observe that, with the current HARPS precision, the 
inclination will be constrained around 2025 if id is close to 10°, 
but if this planet is also coplanar with the other two (/</ ~ 50°), 
its precise value cannot be obtained before 2050. 

than profitable. Indeed, the internal errors of these three instru- 
ments are much higher than those from Keck and HARPS se- 
ries, and they are unable to help in constraining the inclina- 
tions. Moreover, all the radial velocities are relative in nature, 
and therefore each data set included requires the addition of a 
free offset parameter in the orbit fitting procedure, V 'instrument- 

Nevertheless, to be sure that there is no gain in including 
these additional 102 measurements, we performed an adjustment 
to the data using the five instruments simultaneously. The orbital 



parameters corresponding to the best-fit solution are listed in 
Table[3] As expected, this fit yields an adjustment of yjx* = 1.46 
and rms = 3.01 ms _1 , which does not represent an improvement 
with respect to the fit listed in Table|2] The inclination of planet 
c decreases by 2.5°, but this difference is within the 3<x uncer- 
tainty of the best-fit values. Therefore, the orbital parameters de- 
termined only with data from Keck and HARPS will be adopted 
(TableEJ. 

4. Dynamical analysis 

We now analyze the dynamics and stability of the planetary sys- 
tem given in Table|2] Because of the two outermost planets' 
proximity and high values of their masses, both planets are af- 
fected by strong planetary perturbations from each other. As a 
consequence, unless a resonant mechanism is present to avoid 
close encounters, the system cannot be stable. 

4.1. Secular coupling 

As usual in planetary s ystems, there is a strong coupling within 
the secular system (see iLaskanl 1 990h . Therefore, both planets b 
and c precess with the same precession frequency g2, which is 
retrograde with a period of 8.74 years. The two periastron are 
thus locked and Atzr = m c - T&b oscillate around 0° with an am- 
plitude of abou t 25° (Fig .|6l). This b ehavior, a l ready mentioned in 
earlier studies dLaughlin &"C hambers 2001; Lee & Peale 2002; 
Beau ge et alJ 12003). is not a dynamical resonance, but merely 
the result of the linear secular coupling. To present the solution 
in a clearer way, it is then useful to make a linear change o f vari- 
ables into eccentricity and inclination proper modes (see Laskar 
1990). In the present case, the linear transformation is numeri- 
cally obtained by a frequency analysis of the solutions. For the 
dynamical analysis, to understand the evolution of the inclina- 
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tions in this system, we have expressed all the coordinates in the 
reference frame of the invariant plane, orthogonal to the total 
angular momentum of the system. The longitude of node O and 
inclination / of this invariant plane are 

Q = -0.445256°; 7 = 48.767266°. (1) 
Using the classical complex notation, 

z k = e k e^- & = sin(i k /2)e int ; (2) 
for k = b,c, d, we have for the linear Laplace-Lagrange solution 



'z d ) 
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0.139690 0.000843 -0.001265 
-0.000069 0.265688 -0.006501 
0.000063 0.037609 0.006970 
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(3) 



l£d\_l 0.009449 +0.002723 \ I vi 
\£ c J ~ \ -0.000004 -0.013783 J \ v 2 

and the additional approximate relation 
& ~ -0.000028 vj + 0.003301 v 2 . 



(4) 



(5) 



The proper modes in inclination are then given to a good approx- 
imation as \'k « e l(Skt+, '' k -', where su and i/^ are given in Table[4] 

With Eq.[3] it is then easy to understand the meaning of the 
observed libration between the periastrons mi, and tjt c . Indeed, 
for both planets b and c, the dominant term is the u 2 term with 
frequency g 2 . They thus both precess with an average value of 
g 2 . In the same way, both nodes Q/, and Q c precess with the 
same frequency 52- It should also be noted that Eqs.[3]|4] and|5] 
provide good approximations of the long-term evolution of the 
eccentricities and inclinations. Indeed, in Fig. [5] we plot the ec- 
centricity and the inclination with respect to the invariant plane 
of planets b, c, d, with initial conditions from Table [2] At the 
same time, we plot with solid black lines the evolution of the 
same elements given by the above secular, linear approximation. 

The eccentricity and inclination variations are very limited 
and described well by the secular approximation. The eccentric- 
ity of planets b and c are within the ranges 0.028 < < 0.050 
and 0.258 < e c < 0.279, respectively. These variations are 
driven mostly by the rapid secular frequency #3 - g 2 , of period 
2^/(^3 - £2) ~ 8.62 yr (Table|4]i. The eccentricity of planet d is 
nearly constant with 0.136 < ej < 0.142. 

The inclinations of b and c with respect to the invariant plane 
are very small with 0.36° < i b < 0.39° and 1.54° < i c < 1.61°, 
respectively. This small variation is mostly caused by the non- 
linear secular term 2(s2 - £2) of period 4.791 yr. Although the 
inclination of d is not well constrained, using the initial con- 
ditions of Table |2l one finds small variations in its inclination 
0.76° < i c < 1 .42°, which are driven by the secular frequency 
S\ - 52 of period 138.3 yr (Table©. 

With the present 1 1 years of observations covered by Keck 
and HARPS, the most important features that allow us to con- 
strain the parameters of the system are those related to the rapid 
secular frequencies g 2 and s 2 , of periods 8.7 yr and 99 yr, which 
are the precession frequencies of the periastrons and nodes of 
planets b and c. 




10 20 30 40 50 60 
time (yr) 



80 90 100 



CT> 
CD 



where the proper modes are obtained from the Zk by inverting 
the above linear relation. To good approximation, we then have 
Uk ~ e l(gkt+ * k) , where gk and <pk are given in Table[4] For the in- 2- 
clinations, due to the conservation of angular momentum, there c 
are only two proper modes, vi, v 2 , which we can define with the 
invertible linear relation 




100 200 300 400 500 600 700 800 900 1000 
time (yr) 

Fig. 5. Evolution of the GJ 876 eccentricities (top) and inclina- 
tions (bottom) with time, starting with the orbital solution from 
Table|2] The color lines are the complete solutions for the vari- 
ous planets (b: blue, c: green, d: red), while the black curves are 
the associated values obtained with the linear, secular model. 

Table 4. Fundamental frequencies and phases for the orbital so- 
lution in Table|2] 





Frequency 


Period 




Phase 




deg/yr 


yr 




deg 


n b 


2154.322532 


0.167106 




36.0925 


n c 


4349.841605 


0.082762 


4o 


158.1879 


"«/ 


67852.886097 


0.005306 




29.6886 


Si 


1.003492 


358.747259 


fa 


170.0005 


g2 


-41.196550 


8.738596 


<p2 


-86.0030 


g3 


0.555047 


648.593233 


<t>3 


89.8880 


Si 


-1.021110 


352.557387 


<A, 


3.9446 


$2 


-3.624680 


99.319103 




64.1066 


h 


245.918001 


1.463903 


60 


-48.6812 



n b , n c and n d are the mean motions, g\, g% and g 3 the main secular 
frequencies of the periastrons, s\ and si the main secular frequencies of 
the nodes, and l e the libration frequency of the resonant angles. 



4.2. The 2:1 mean motion resonance 

The ratio of the orbital periods of the two outermost planets de- 
termined by the fitting process (Tabled is Ph/P c = 2.018, sug- 
gesting that the system may be trapped in a 2:1 mean motion 
resonance. This resonant motion has already been reported in 
previous works (jLaughlin & Chambersl200ltlRivera & Lissauerl 
120011: iLee & Pealell2002t iRivera et al.1 120051) . To test the accu- 
racy of this scenario, we performed a frequency analysis of 
the orbital solution listed in Table|2] computed over 100 kyr. 
The orbits of the three p lanets are integrated with the symplec- 
tic integrator SABA4 of lLaskar & Robuteil d2001h . using a step 
size of 2 x 10~ 4 years. We conclude that in the nominal solu- 
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Fig. 6. Variation in the resonant argument 6 = 2 At, -A c -m* 2 with 
time (top). The difference Aa> = mt - m c also oscillates around 
0° with a 25° amplitude. 

tion of Table[2] planets b and c in the GJ 876 system indeed 
show a 2:1 mean motion resonance, with resonant arguments: 
0b — 2Ab — A c — ETj and G c — 2Aj — A c — m c . 

If we analyze these arguments, it is indeed difficult to dis- 
entangle the proper libration from the secular oscillation of the 
periastrons angles uJb,Tn c . It is thus much clearer to switch to 
proper modes Uk, V* by means of the linear transformation ([3]|4]l. 
Indeed, if = & the argument 

6 = 2A b -A c + td\ (6) 

is in libration around 0° with a very small total amplitude of 3.5° 
(Pig.[6j» with only 1.8° amplitude related to the proper libration 
argument 9 with period 2n/lg = 1.46 yr. The remaining terms 
of Table[5]are forced terms related to planetary interactions. The 
fundamental frequencies associated with this libration argument 
(Table© indeed verify the resonance relation 2nb-n c -g2 — up 
to the precision of the determination of the frequencies (» 1CT 5 ). 

4.3. Coplanar motion 

The orbits of the planets in the Solar System are nearly coplanar, 
that is, their orbital planes remain within a few degrees of incli- 
nation from an inertial plane perpendicular to the total angular 
momentum of the system. The highest inclination is obtained for 
Mercury (/ ~ 7°), which is also the smallest of the eight planets. 

A general belief that planetary systems would tend to resem- 
ble the Solar System persists, but there is no particular reason for 
all planetary systems to be coplanar such as our own. Star forma- 
tion theories require an accretion disk in which planets form, but 
close encounters during the forma tion process ca n increase the 
eccentricities and inclinations (e.g. lLee et alj2007h . That planets 
are found to have significant eccentricity values is also consistent 



Table 5. Quasi-periodic decomposition of the resonant angle 
8 = 2Ab — Ac — tu\ for an integration over 100 kyr of the orbital 
solution in Table|2| 



Combination v ; A, <p t 



n b 


n c 


82 


£3 


k 


(deg/yr) 


(dog) 


(cleg) 














1 


245.9180 


1.810 


-48.681 


1 


1 


-2 








6586.5572 


0.567 


96.286 





1 


-1 








2195.5191 


0.569 


32.095 








-1 


1 





41.7516 


0.411 


-94.109 


1 





-1 








4391.0382 


0.255 


-25.809 


2 


1 


-3 








10977.5954 


0.156 


-19.523 








1 


-1 


1 


204.1664 


0.120 


-44.572 



We have 8 = Aj cos(y, t + <pj). We only provide the first 7 largest 
terms that are identified as integer combinations of the fundamental fre- 
quencies given in Table|4] 




50 100 150 200 

time (year) 

Fig. 7. Evolution of the mutual inclination of the planets b and 
c. Since the maximal inclination between the two planets is 2.0°, 
we conclude that the orbits are almost coplanar, as in our Solar 
System. 



with their having high inclinations (e.g. lLibert & Tsiganisl2 009). 
In particular, studies of Kuiper belt objects indicate that there 
is an important inclinati on excitation when the bodies sweep 
secular resonances (e.g. [Nagasaw a et al.l l2000t iLi et al . 2008), 
a mechanism that also appe ars to be applicable to planets (e.g. 
iThommes & Lissauerll2003l) . 

The question of whether extra-solar planetary systems 
are also nearly coplanar or not is thus important. Although 
about 45 extra-solar multi-planet systems have been reported, 
their true inclinations have so far been determined in only 
two cases. One is the planetary system around the pulsar 
PSR B 1257 +12, discovered b y precise timing measurements of 
pulses (IWolszczan & Frail|l992|). for which t he orbits seem to be 
almost coplanar dKonacki & Wolszczanll2003l) . The other case is 
GJ 876, the only planetary system around a main-sequence star 
for which one can access the inclinations. This result is possi- 
ble because of the large amount of data already available and 
because the two giant planets show strong gravitational interac- 
tions. 

The analytical expressions and the plots versus time of the 
orbital elements in reference frame of the invariant plane show 
that the inclinations of planets b and c are both very small 
(Eqs.[4][5] Fig.[5j. To test the coplanarity of this system, we plot 
in Fig. [7] the evolution of the mutual inclination with time. We 
verify that / = 1.958° + 0.045°, and hence conclude that the 
orbits of planets b and c are indeed nearly coplanar. 
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In our best-fit solution (Tabled, we also assumed that the 
orbit of planet d is nearly coplanar with the other two planets, 
but there is no clear physical reason that justifies our choice. 
Indeed, its orbit should be more perturbed by the oblateness of 
the st ar, rather than the other planets (Goldreich 19651 ICorreial 
l2009h . It is therefore more likely that planet d orbits in the equa- 
torial plane of the star. This plane can be identical to the orbital 
plane of the outer planets or not, depending on the evolutionary 
process. According to Fig.[2l with the current HARPS precision 
of ~1 m/s, we will have to wait until about 2050 to confirm an 
inclination close to 50°. 

Adopting Qd = 0°, and initial values of id ranging from 10° 
to 90°, we found best-fit solutions with very similar values of 
reduced yfx* as the solution in Table [2] For inclinations lower 

than 10°, the best-fit solution becomes worse {^Jx 2 = 1.58 for 
id = 5°) and infers low eccentricity values for planet d. In addi- 
tion, planets b and c are no longer nearly coplanar. We therefore 
conclude that a lower limit to the inclination of planet d can be 
set to be around 10°. On the other hand, a lack of transit detec- 
tions only allows inclina tions close to 90° if the planet is very 
dense dRivera et al.ll2005l) . 

4.4. Stability analysis 

To analyze the stability of the nominal solution (Table|2| and 
confirm the presence of t he 2:1 reson ance, we performed a 
global frequency analysis (Laskar Hl993l) in the vicinity of this 
solution (Fig. lHi, in the same w ay as achieved for the HD 45364 
system bv lCorreia eTaf] d2009h . 

For each planet, the system is integrated on a regular 2D 
mesh of initial conditions, with varying semi-major axis and ec- 
centricity, while the other parameters are retained at their nom- 
inal values (Tabled. The solution is integrated over 800 yr for 
each initial condition and a stability indicator is computed to be 
the variation in the measured mean motion over the two con- 
secutive 400 yr intervals of time. For regular motion, there is 
no significant variation in the mean motion along the trajectory, 
while it can vary significantly for chaotic trajectories. The re- 
sult is reported in color in Fig. [8] where "red" represents the 
strongly chaotic trajectories, and "dark blue" the extremely sta- 
ble ones. To decrease the computation time, we averaged the 
orbit of planet d oy er its mean-motion and periastron, following 
iFarago et al.ld2009l) . 

In the two top plots (Figs. |8^,b), the only stable zone that 
exists in the vicinity of the nominal solutio n corresponds to the 
stable 2:1 reso nant areas. As for HP 45 364 (ICorreia et al.ll2009h 
and HD 60532 (lLaskar & Correial2 009) planetary systems, there 
is a perfect coincidence between the stable 2:1 resonant islands, 
and curves of minimal y 1 obtained by comparison with the ob- 
servations. Since these islands are the only stable zones in the 
vicinity, this picture presents a very coherent view of dynamical 
analysis and radial velocity measurements, which reinforces the 
confidence that the present system is in a 2: 1 resonant state. 

The scale of the two bottom panels shows the remarkable 
precision of the data in light of the dynamical environment of 
the system. The darker structures in these plots can be identified 
as secondary resonances. For instance, in the bottom left plot 
of Fig. 7, the two dark horizontal lines starting at e c =2 0.269 
and 0.260 correspond, respectively, to lg + 6g2 - 6g3 + 51-52 
and lg + 6g2 - 2gi - 4s2- We can also barely see the secondary 
resonance lg + 6g2 - 6gi just above the former one. In the bottom 
right plot, the bluer line starting at e\, 0.023 corresponds to lg + 
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Fig. 8. Stability analysis of the nominal fit (Tabled of the 
GJ 876 planetary system. For a fixed initial condition of planet 
b (left) and planet c (right), the phase space of the system is 
explored by varying the semi-major axis a* and eccentricity ej 
of the other planet, respectively. The step size is 0.0002 AU in 
semi-major axis and 0.005 in eccentricity in the top panels and 
2x 10~ 6 AU and 2xl0~ 4 in bottom panels. For each initial condi- 
tion, the syste m is integrated ove r 800 yr with an averaged model 
for planet d (Farag o et al.ll2009h and a stability criterio n is de- 
rived with th e freq uency analysis of the mea n longitude (Laskar 
ll99QLl!993h . As in lCorreia et all d2005l I2009I the chaotic diffu- 
sion is measured by the variation in the frequencies. The "red" 
zone corresponds to highly unstable orbits, while the "dark blue" 
region can be assumed to be stable on a billion-years timescale. 
The contour curves indicate the value of x 2 obtained for each 
choice of parameters. It is remarkable that in the present fit, there 
is an exact correspondence between the zone of minimal x 1 and 
the 2: 1 stable resonant zone, in "dark blue". 



6g2 - 4g3 - 252. A longer integration time would most probably 
provide more details about the resonant and secular dynamics. 

4.5. Long-term orbital evolution 

From the previous stability analysis, it is clear that planets b and 
c listed in Table|2]are trapped in a 2:1 mean motion resonance 
and stable over a Gyr timescale. Nevertheless, we also tested 
this by performing a numerical integration of the orbits using 
the symplectic integrator SABA4. 

Because the orbital period of the innermost planet is shorter 
than 2 day, we performed two kinds of experiments. In the 
first one, we directly integrated the full planetary system over 
100 Myr with a step size of 2 x 10~ 4 years. Although tidal ef- 
fects may play an important role in this system evolution, we 
did not include them. The result is displayed in Fig. [9] showing 
that the orbits indeed evolve in a regular way, and remain stable 
throughout the simulation. For longer timescales, we needed to 
averag e the orbit of the inner planet as explained in Farag oet al.l 
(2009). We then use a longer step size of 2x 10 3 years and inte- 
grated the system over 5 Gyr, which corresponds to the maximal 
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Fig. 10. Possible location of an additional fourth planet in the GJ 876 system. The stability of an Earth-size planet (K — 1 m/s) is 
analyzed, for various semi-major axis versus eccentricity (top), or mean anomaly (bottom). All the angles of the putative planets 
are set to 0° (except the mean anomaly in the bottom panels), its inclination to 50°, and in the bottom panels, its eccentricity to 0. 
The stable zones where additional planets can be found are the dark blue regions. Between the already known planets, we can find 
stability in a small zone corresponding to a 4: 1 mean motion resonance with planet b (and 2: 1 resonance with planet c). White lines 
represent the collisions with the already existing planets, given by a cross. 




-0.2 -0.1 0.1 0.2 

x(AU) 



Fig. 9. Long-term evolution of the GJ 876 planetary system over 
100 Myr starting with the orbital solution from Table|2] We did 
not include tidal effects in this simulation. The panel shows a 
face-on view of the system invariant plane, x and y are spatial 
coordinates in a frame centered on the star. Present orbital solu- 
tions are traced with solid lines and each dot corresponds to the 
position of the planet every 5 kyr. The semi-major axes (in AU) 
are almost constant (0.209 < a h < 0.214; 0.131 < a c < 0.133 
and 0.02110 < a d < 0.02111), and the eccentricities present 
slight variations (0.028 < e b < 0.050; 0.258 < e c < 0219 and 
0.136 < e d < 0.142). 



estimated age of the central star. The subsystem consisting of the 
giant planets b and c remained stable. 

In spite of the strong gravitational interactions between the 
two planets locked in the 2:1 mean motion resonance, both or- 
bital eccentricities and inclinations exhibit small variations that 
are mostly driven by the regular linear secular terms (Fig|5j- 
These variations occur far more rapidly than in our Solar System, 
which enabled their direct detection using only 11 yr of data 
taken with Keck and HARPS together. It is important to notice, 
however, that for initial inclinations of planet d very different 
from 50°, the orbital perturbations from the outer planets may 
become important. For instance, using an orbital solution with 
id = 90°, we derive an eccentricity and inclination variations of 
0.10 <e d < 0.35 and 10° < i d < 90°. 

4.6. Additional constraints 

The stability analysis summarized in Fig. [8] shows good agree- 
ment between the 2: 1 resonant islands and the^ 2 contour curves. 
We can thus assume that the dynamics of the two known plan- 
ets is not disturbed much by the presence of an additional planet 
close-by. The same is true for the innermost planet, which has 
an orbital period of 2 day, since the gravitational interaction with 
the parent star is too strong to be destabilized. 

We then tested the possibility of an additional fourth planet 
in the system by varying the semi-major axis, the eccentricity, 
and the longitude of the periastron over a wide range, and per- 
forming a stability analysis (Fig. ITOb . The test was completed 
for a fixed value K — 1 m/s, corresponding to an Earth-size ob- 
ject. We also performed a simulation of a Neptune-size object 
(K = 10 m/s) without significant changes in its dynamics. In this 
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last case, however, an object of this size would have already been 
detected in the data. 

From this analysis (Fig. [TOl . one can see that stable orbits 
are possible beyond 1 AU and, very interestingly, also for orbital 
periods around 15 days (~0.083 AU), which correspond to bod- 
ies trapped in a 4: 1 mean motion resonance with planet b (and 
2: 1 resonance with planet c). Among the already known planets, 
this is the only zone where additional planetary mass compan- 
ions can survive. With the current HARPS precision of ~1 m/s, 
we estimate that any object with a minimum mass M > 2M @ 
would already be visible in the data. Since this does not seem to 
be the case, if we assume that a planet exists in this resonant sta- 
ble zone, it should be an object smaller or not much larger than 
the Earth. 

The presence of this fourth planet, not only fills an empty gap 
in the system, but can also help us to explain the anomalous high 
eccentricity of planet d (ej ~ 0.14). Indeed, tidal interactions 
with the star sho uld have circular ized the orbit in less than 1 Myr, 
but according to Mardling (2007), the presence of an Earth-sized 
outer planet may delay the tidal damping of the eccentricity. 

We propose that additional observational efforts should be 
made to search for this planet. The orbital period of only 15 day 
does not require long time of telescope and a planet with 
the mass of the Earth is at reach with the present resolution. 
Moreover, assuming coplanar motion with the outermost plan- 
ets, we obtain the exact mass for this planet and not a mini- 
mal estimation. We have reanalyzed the residuals of the system 
(Fig-Ell, but so far no signal around 15 days appears to be present 
with the current precision. However, since this planet will be in 
a 2:1 mean motion resonance with planet c, we cannot exclude 
that the signal is hidden in t he eccentricity of the giant planet 
dAnglada-Escude et alj|2008h . 

5. Discussion and conclusion 

We have reanalyzed the planetary system orbiting the star 
GJ 876, using the new high-precision observational data ac- 
quired with HARPS. Independently from previous observations 
with other instruments, we can confirm the presence of planet d, 
orbiting at 1.93785 days, but in an elliptical orbit with an eccen- 
tricity of 0.14 and a minimum mass of 6.3 M s . 

By combining the HARPS data with the data previously 
taken at the Keck Observatory, we are able to fully characterize 
the b and c planets. We find i b =48.9°±1.0° andi c = 48.1°±2.1°, 
which infers for the true masses of the planets Mb = 2.64 Mj up 
and M c = 0.83 Mj up , respectively. We hence conclude that the 
orbits of these two planets are nearly coplanar. The gravitational 
interactions between the outer planets and the innermost planet 
d may also allow us to determine its orbit more accurately in 
the near future. With the current precision of HARPS of ~1 m/s 
for GJ 876, we expect to detect the true inclination and mass of 
planet d within some decades. 

A dynamical analysis of this planetary system confirms that 
planets b and c are locked in a 2: 1 mean motion resonance, which 
ensures stability over 5 Gyr. In the nominal solution, the reso- 
nant angles 6 b = 2Ab - A c + Wb and 9 C = 2Ab - A c + m c are in 
libration around 0°, which means that their periastrons are also 
aligned. This orbital configuration may have been reached by 
means of the dissipative process of plane t migration durin g the 
early stages of the system evolution (e.g. Crida et al. 2008). By 
analyzing the proper modes, we are able to see that the ampli- 
tude libration of the proper resonant argument 9 - 2Ab - A c + vj* 2 
can be as small as 3.5°, of which only 1.8° is related to the libra- 
tion frequency. It is thus remarkable that the libration amplitude 



is so small, owing to the width of the stable resonant island be- 
ing large. We can thus assume that the libration amplitude has 
been damped by some dissipative process. This singular plane- 
tary system may then provide important constraints on planetary 
formation and migration scenarios. 

Finally, we have found that stability is possible for an Earth- 
size mass planet or smaller in an orbit around 15 days, which is 
in a 4: 1 mean motion resonance with planet b. The presence of 
this fourth planet, not only fills an empty gap in the system, but 
can also help us to explain the anomalous high eccentricity of 
planet d (ej ~ 0.14), which should have been damped to zero 
by tides. Because of the proximity and low mass of the star, a 
planet with the mass of the Earth should be detectable at the 
present HARPS resolution. 

We conclude that the radial-velocity technique is self suffi- 
cient for fully characterizing and determining all the orbital pa- 
rameters of a multi-planet system, without needing to use as- 
trometry or transits. 
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